home *** CD-ROM | disk | FTP | other *** search
/ Precision Software Appli…tions Silver Collection 4 / Precision Software Applications Silver Collection Volume 4 (1993).iso / stats / chadyn.exe / YLORENZ.C < prev    next >
C/C++ Source or Header  |  1988-11-28  |  7KB  |  272 lines
  1. /********************************* YLorenz.C *********************************/
  2. /*********************** (C) 1986,7,8 by JAMES A. YORKE **********************/
  3.  
  4. #include "yinclud.h"
  5.  
  6.  
  7. /*********** Lorenz  1 mode --- 3 equations  ***********/
  8. DELorenz(k, Y)
  9. double  Y[],
  10.         k[];
  11. {
  12.     int     i,
  13.             base;
  14.     double  x1,
  15.             y1val,        /* Y1 is a global variable */
  16.             z1;
  17.     double  J00,
  18.             J01,
  19.             J10,
  20.             J11,
  21.             J12,
  22.             J20,
  23.             J21,
  24.             J22;
  25.     x1 = Y[0];
  26.     y1val = Y[1];
  27.     z1 = Y[2];
  28.  
  29.     k[0] = -sigma * (x1 - y1val);
  30.     k[1] = rho * x1 - y1val - x1 * z1;
  31.     k[2] = x1 * y1val - beta * z1;
  32.     k[3] = 1;        /* time */
  33.     if(num_lyap > 0) {
  34.         J00 = -sigma;
  35.         J01 = sigma;
  36.         J10 = rho - z1;
  37.         J11 = -1.;
  38.         J12 = -x1;
  39.         J20 = y1val;
  40.         J21 = x1;
  41.         J22 = -beta;
  42.     }
  43.     for(i = 0; i < num_lyap; i++) {
  44.         base = lyapzero + vec_dim * i;
  45.         k[base + 0] =
  46.             J00 * Y[base] + J01 * Y[base + 1];
  47.         k[base + 1] =
  48.             J10 * Y[base] + J11 * Y[base + 1] + J12 * Y[base + 2];
  49.         k[base + 2] =
  50.             J20 * Y[base] + J21 * Y[base + 1] + J22 * Y[base + 2];
  51.     }
  52. }
  53.  
  54.  
  55. initLorenz() {
  56.  /* Lorenz */
  57.     int     DELorenz();
  58.     vec_dim = 3;        /* the dimension of the Lyapunov vectors =
  59.                    phase space dim */
  60.     num_lyap = 0;        /* the number of Lyapunov numbers to be
  61.                    computed <= vec_dim  */
  62.     lyapzero = 4;        /* y[lyapzero] is the zeroth coord of the
  63.                    zeroth lyapunov vector */
  64.     dim = lyapzero + num_lyap * vec_dim;
  65.  /* needed for rungekutta */
  66.  
  67.     X_upper = 25.0;        /* x scale */
  68.     X_lower = -25.;
  69.     Y_upper = 60.;        /* y scale */
  70.     Y_lower = 0.;
  71.  
  72.     Y_coord = 2;
  73.  
  74.     sigma = 10;
  75.     rho = 28.0;
  76.     beta = 8./ 3.;
  77.  
  78.     step = 1./ 100.;
  79.     DE = DELorenz;        /* DE is a pointer to a function   */
  80. }
  81. DESam(k, Y)
  82. double  Y[],
  83.         k[];
  84. {
  85.     int     i,
  86.             base;
  87.     double  x1,
  88.             y1val,        /* Y1 is a global variable */
  89.             z1;
  90.     double  J00,
  91.             J01,
  92.             J02,
  93.             J10,
  94.             J11,
  95.             J12,
  96.             J20,
  97.             J21,
  98.             J22;
  99.     x1 = Y[0];
  100.     y1val = Y[1];
  101.     z1 = Y[2];
  102.  
  103.     k[0] = C1 * x1 + C2 * y1val - z1;
  104.     k[1] = C3 * x1 + C5 * y1val - C6 * x1 * x1 * x1;
  105.     k[2] = rho * x1;
  106.     k[3] = 1;        /* time */
  107.     if(num_lyap > 0) {
  108.         J00 = C1;
  109.         J01 = C2;
  110.         J02 = -1;
  111.         J10 = C3 - 3 * C6 * x1 * x1;
  112.         J11 = C5;
  113.         J12 = 0;
  114.         J20 = rho;
  115.         J21 = 0;
  116.         J22 = 0;
  117.     }
  118.     for(i = 0; i < num_lyap; i++) {
  119.         base = lyapzero + vec_dim * i;
  120.         k[base + 0] =
  121.             J00 * Y[base] + J01 * Y[base + 1] + J02 * Y[base + 2];
  122.         k[base + 1] =
  123.             J10 * Y[base] + J11 * Y[base + 1] + J12 * Y[base + 2];
  124.         k[base + 2] =
  125.             J20 * Y[base] + J21 * Y[base + 1] + J22 * Y[base + 2];
  126.     }
  127. }
  128.  
  129.  
  130. initSam() {
  131. #ifdef IGNORE
  132.     TESTMAP2 ("sam")
  133.         fputs(
  134.             "SAM  Samardzia system: x'= C1*x+C2*y-z, y'= C3*x+C5*y-C6*x^3, z'= rho*x\n"
  135.             ,output);
  136.     TESTMAP("sam")
  137.         fputs(
  138.             "Numbering of variables: y[0]=x ,y[1]=y, y[2]=z, y[3]=time                  \n"
  139.             ,output);
  140. #endif
  141.  /* Nick Samardzia - DuPont 302-695-3009 */
  142.  
  143.     int     DESam();
  144.  
  145.     num_lyap = 0;        /* the number of Lyapunov numbers to be
  146.                    computed <= vec_dim  */
  147.     vec_dim = 3;        /* the dimension of the Lyapunov vectors =
  148.                    phase space dim */
  149.     lyapzero = 4;        /* y[lyapzero] is the zeroth coord of the
  150.                    zeroth lyapunov vector */
  151.     X_upper = 2.0;        /* x scale */
  152.     X_lower = -2.;
  153.     Y_lower = -2.;        /* y scale */
  154.     Y_upper = 2.;
  155.     X_coord = 0;
  156.     Y_coord = 1;
  157.  
  158.     C1 =.032;
  159.     C2 =.5;
  160.     C3 = 1.;
  161.     C5 = -.1;
  162.     C6 = 1.;
  163.     rho = 0.1;
  164.     step =.04;
  165.     rho_step =.02;
  166.     DE = DESam;        /* DE is a pointer to a function   */
  167. }
  168.  
  169.  
  170. ReturnMap() {            /* this routine will iterate (*DEsolver) until
  171.                    Fn() changes sign from positive to negative,
  172.                    that is, until old_Fn>0 and new_Fn<=0); the
  173.                    time between crossings(= "time" at the end
  174.                    of this routine) is returned; for
  175.                    differential equations it is necessary to
  176.                    check to see if some variable(possibly
  177.                    time) should be taken mod something; Fn()
  178.                    must be chosen so that the mod operation
  179.                    does not decrease Fn from positive to
  180.                    negative; if the cross section is going to
  181.                    be set up so that time t = some t0, and t is
  182.                    increasing and is taken mod 1, say, then
  183.                    define Fn() = t0 - t so that the mod
  184.                    increases Fn; intermediate data on finding
  185.                    cross section is printed out when "printer"
  186.                    == 3 */
  187.     double  time,
  188.             savestep,
  189.             old_Fn,
  190.             new_Fn,
  191.             old_y[EQUATIONS];
  192.  
  193.     new_Fn = (*Fn) ();
  194.     old_Fn = new_Fn;    /* since new_fn = old_fn, the "for()" will be
  195.                    traversed at least once */
  196.     for(time = 0; (level >= PROCESS) && ((old_Fn <= 0) || (new_Fn > 0));
  197.             time = time + step) {
  198.                 /* continue cycling unless Fn() decreases from
  199.                    > 0 to <= 0 */
  200.         old_Fn = new_Fn;
  201.         store(old_y, y, dim);
  202.         (*DEsolver) ();
  203.         (*modPointer) ();
  204.                 /* this = null() for non differential equations
  205.                    and for many differential equations */
  206.     /* if this mod causes(*Fn)() to decrease past 0, we have a problem */
  207.         new_Fn = (*Fn) ();
  208.         if(printer == 3) {
  209.             scr_rowcol(2, 0);
  210.             printf(
  211.                     "quick: Hit 2 to stop all this printing; then hit c to clear the screen   \n");
  212.             Interrupt();
  213.             printf(
  214.                     "(*Fn)() = %12.10lf    time = %12.10lf, step = %12.10lf         \n"
  215.                     ,new_Fn, time, step);
  216.         }
  217.     }
  218.  
  219.     savestep = step;
  220.     while(level >= PROCESS && fabs (new_Fn) >.00000001) {
  221.                 /* fabs() is absolute value of a double
  222.                    precision number */
  223.         step = step * old_Fn / (old_Fn - new_Fn);
  224.                 /* estimated time step needed to hit(*Fn)() =
  225.                    0 starting from old_y  */
  226.         sixth_step = step / 6;/* for rungekutta() */
  227.         halfstep = step / 2;
  228.  
  229.  
  230.         store(y, old_y, dim);
  231.         (*DEsolver) ();
  232.         (*modPointer) ();
  233.                 /* this = null() for non differential equations
  234.                    and for many differential equations */
  235.     /* if this modulo causes(*Fn)() to decrease past 0, we have a problem 
  236.     */
  237.         new_Fn = (*Fn) ();
  238.         if(printer == 3) {
  239.             scr_rowcol(2, 0);
  240.             printf(
  241.                     "REFINING CROSSING POINT: new Fn = %12.10lf; new step = %12.10lf  \n"
  242.                     ,new_Fn, step);
  243.             Interrupt();
  244.         }
  245.     }
  246.     time = time - savestep + step;
  247.                 /* trajectory overshot in first part so what
  248.                    was then step must be subtracted off and we
  249.                    add the partial step which at this moment is
  250.                    called step */
  251.     step = savestep;
  252.     sixth_step = step / 6;    /* for rungekutta() */
  253.     halfstep = step / 2;
  254.     if(printer == 3)
  255.         printf("time for loop: %12.10lf     \n", time);
  256.     return;
  257. }
  258.  
  259. initLRetMap() {
  260.     extern double   Fn3 ();
  261.     extern int      ReturnMap();
  262.  
  263.     Fn = Fn3;        /* for poincare return map */
  264.     iteratee = ReturnMap;
  265.     initLorenz();
  266.     Y_coord = 1;        /* X_coord = 0 */
  267.     Y_upper = 10.;        /* x & y scales are symmetric */
  268.     Y_lower = -10.;
  269.     X_upper = 10.0;        /* x scale */
  270.     X_lower = -10.;
  271. }
  272.